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Based on fundamental measure theory, a Helmholtz free energy density functional for three- 
component mixtures of hard spheres with general, non-additive interaction distances is constructed. 
The functional constitutes a generalization of the previously given theory for binary non-additive 
mixtures. The diagrammatic structure of the spatial integrals in both functionals is of star-like (or 
tree-like) topology. The ternary diagrams possess a higher degree of complexity than the binary 
diagrams. Results for partial pair correlation functions, obtained via the Ornstein-Zernikc route 
from the second functional derivatives of the excess free energy functional, agree well with Monte 
Carlo simulation data. 
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I. INTRODUCTION 

Rosenfeld's fundamental measures theory (FMT) for 
additive hard sphere mixtures 0, has become a cor- 
nerstone of classical density functional theory (DFT) Q . 
FMT features in numerous applications to a wide vari- 
ety of interesting phenomena in liquids ■ I n additive 
hard sphere mixtures the interaction distance between 
unlike components is taken to be the arithmetic mean 
of the (like-species) diameters. Additive mixtures are of- 
ten considered as prototypical in the description of liquid 
mixtures. Nevertheless, non-additivity is a generic fea- 
ture p} that arises very naturally in effective interactions, 
e.g. due to the depletion effect Q, or when integrating out 
solvent degrees of freedom in electrolytes [H, |T(| ■ A gen- 
eralization of FMT to binary non-additive hard sphere 
(NAHS) mixtures was based on the scalar version 
of the additive hard sphere functional Q, and has been 
used successfully in the investigation of bulk [l2rll3 | and 
of interfacial [1 51 ] phenomena that occur in NAHS mix- 
tures. 

Recently significant progress has been made in the for- 
malization of the mathematical structure of FMT. This 
includes i) insights into the geometry of the non-local 
aspects of the theory, i.e. the algebraic group structure 
and symmet ry prop erties of the convolution kernel ma- 
trix K(R, r) [llj, lla [l7| , which controls the range of non- 
locality in the functional. Here R is a fixed lengthscale 
and r is the radial distance in three-dimensional space. 
Formally, K(i?, r) is a 4 x 4-matrix that is indexed by 
powers of lengthscale. Its remarkable algebraic proper- 
ties [IH allow to view it as an object that is suitable to 
add or remove a layer of thickness R from a given sphere. 
Furthermore, ii) the FMT for additive hard sphere mix- 
tures was obtained from a tensorial-diagrammatic series 
in density [18J. The diagrams possess star-like topology, 
which is an approximation of the combinatorial complex- 
ity of the exact virial expansion. Apart from a central 
space integral, all field points are integrations over the 
density field(s), as they are in the exact virial expansion. 



The bonds, however, are weight functions rather than 
Mayer functions, and possess only half the range of the 
(hard sphere) Mayer function. In one dimension, the re- 
sult from the series is equal to Percus' exact functional 
[lj|. In three dimensions it gives the Kierlik-Rosinberg 
form of FMT. The five-dimensional hard hypersphere 
version is investigated in detail in Ref. (l8| . Comparison 
to data from the literature for bulk structure and ther- 
modynamics demonstrates the capability of this theory 
for study bulk and inhomogeneous hypersphere mixtures. 



In the present work we generalize the tensorial- 
diagrammatic series to non-additive hard sphere inter- 
actions. We use the kernel matrix K(i?,r) as a further 
type of bond. This enables us to represent the binary 
NAHS functional of Ref. as a series of diagrams that 
are formed by two stars, one for each species. Here the 
center of one star is connected to the center of the sec- 
ond star by a K-bond. Formulating an FMT for gen- 
eral ternary non-additive hard sphere mixtures requires 
to modify this topology. Treating special cases of ternary 
mixtures that have a suitably high degree of additi vity 
amounts to straightforward generalizations of Ref. [111 ]. 
This is the case when the cross interaction between (say) 
species 1 and 2 is additive, and the non-additivities be- 
tween 13 and 23 are coupled in a certain way (see the 
discussion below (TT])), and hence cannot be chosen in- 
dependently from each other. The general case, how- 
ever requires the introduction of diagrams with different 
topology. Here three stars, one for each species, are con- 
nected to a central three-arm star. As laid out in detail 
below, all (four) inner junctions are bare space integrals, 
without multiplication by a one-body density. Only the 
outer ends carry multiplication by a (bare) density vari- 
able. We show that the theory predicts the bulk struc- 
ture of ternary mixtures, via the Ornstein-Zernike route, 
with good quality as compared to Monte Carlo simula- 
tion data. 
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II. NON-ADDITIVE HARD SPHERE 
INTERACTIONS 



Non-additive hard sphere mixtures possess pair po- 
tentials Vij (r) between species i and j as a function of 
the center-center distance r of the two particles, that 
are given as Vij(r) = oo if r < o~ij, and zero other- 



wise. In the general case all <Xy are independent of 
each other, except for the trivial symmetry aij = Uji. 
Non-additivity parameters are conventionally defined as 
Aij = 2a l . J /{a li + a-jj) - 1 for i ^ j, where A VJ > -1. 



For additive mixtures all A; 



0. For positive non- 



additivity, > 0, the unlike components interact at a 
larger distance than the arithmetic mean of their diam- 
eters. For negative non-additivity, Ajj < 0, the interac- 
tion distance is smaller than the mean of the diameters. 
While a binary NAHS mixture is characterized by either 
positive or negative non-additivity, in ternary systems 
mixed cases are possible, where not all of the Ay have 
the same sign. In ternary NAHS mixtures the equation 
of state [13] was considered and phase stability was inves- 
tigated using integral equation theory [21|. Phase equi- 
librium was also considered in polydisperse non-additive 
hard sphere systems [HI, HH . 



III. DENSITY FUNCTIONAL THEORY 

A. Overview and choice of lengthscales 

In order to construct an FMT for ternary NAHS mix- 
tures, we first identify suitable lengthscales. Consistent 
with Rosenfeld's additive case [l|, we use the particle 
radii Ri = au/2oi each species i = 1,2,3. For the bi- 
nary mixture 11|, the cross diameter between species 1 
and 2 was decomposed as o~\2 = Ri + R12 + -R27 where 
the lengthscale R12 = er 12 — (011 + 022 )/2 accounts for 
the non-additivity. For the ternary system, we generalize 
this to a decomposition of the cross diameters into sums 
of four contributions, 



o~ij — Ri + di + dj + Rj, 



(1) 



where ij — 12,13,23, and the three lengthscales di, 
i = 1,2,3, control the degrees of non-additivity. Invert- 
ing Eq. (TT]) is possible for ternary mixtures and yields 
dj = (uij + <jjk — Oik — <J jj)/2; here ijk — 123 or any 
permutation thereof. As a simple check, a counting exer- 
cise assures us that the number of parameters Ri (three) 
and di (also three) is enough to represent the six inde- 
pendent components of a^. The binary case above is 
recovered if we set R12 = d\ + c?2- Here the relative split- 
ting of R12 into g?i and d2 is arbitrary. If we keep three 
species, and set d\ = d2 = 0, then only ^3 ^ remains 
in order to control the non-additivities between 13 and 
between 23. For mixtures with such restricted degree 
of non-additivity, as mentioned above, the binary NAHS 
functional can be generalized easily. Figure [1] illustrates 



the different types of decomposition of the cry in binary 
additive (a), binary non-additive (b, c), and ternary non- 
additive (d) mixtures. Note that the di can be negative, 
but that R t + di + dj + Rj > must hold due to (TTJ). 
Furthermore certainly Ri > 0. 
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FIG. 1: Illustration of the decomposition {TJ of the hard 
sphere interaction distances aij into like-species particle radii 
Ri and non-additivity distances di. Thin solid lines repre- 
sent distances between a space point (dots) and a junction; 
thick solid lines represent the distance between two junctions. 
The dashed arrows indicate paths that possess length cy (as 
indicated); the sense of direction is only a guide to the eye. 
Decompositions are shown for the binary additive mixture (a), 
the binary non-additive mixture with the cross diameter be- 
ing divided into three (b) and four (c) contributions, and the 
ternary non-additive mixture (d). The two binary cases (b) 
and (c) are equivalent when R12 = d\ + d-z. 

As a means to control the different lengthscales in the 
densit y fu nctional, we use the kernel matrix K(i?,r) of 
Refs. ML [lj| H3]- Two of its properties render this a 
suitable object for the construction of the density func- 
tional: i) the four scalar Kierlik-Rosinberg weight func- 
tions (where R is identified with the particle radius) 
feature as components and ii) matrices can be chained, 
K(R + R',r) = K(R, r) *K(R',r), where the asterisk de- 
notes the three-dimensional spatial convolution and ma- 
trix multiplication is implied on the right hand side. The 
three-dimensional Fourier transform of K(i?,r) can be 
expressed as K(R,q) — exp(i?G), where q is the radial 
distance in reciprocal space, and the (generator) matrix 
G depends on q and is defined by its components Gj = 
TG^ = 87r,G 3 2 = l^ 2 = -q 2 /(Aw),G 3 = ~q 4 /(8n)- 
all other G,, 1 ' = 0. Here the lower index indicates the 
row and the upper index indicates the column; all Greek 
indices run from to 3 here and in the following. Due 
to the symmetry K^(R,r) = K 3 ^_~ V (R, r) [13, the ma- 
trix K(R,r) has ten independent components [H, Il7j . 
The four scalar Kierlik-Rosinberg [H weight functions 
w u (R,r) are contained herein, i.e. K u (R,r) = w v {R,r), 
explicitly given in real space as Wa(R,r) = Q(R — r), 
w 2 (R,r) = S(R-r), un(R,r) = S'(R - r)/{8%), and 
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w (R,r) = -5"{R - r)/(87r) + 5'{R - r)/(2nr). An ex- 
plicit real-space expression for all further components of 
K(_R, r) can be found in fill ]. 

The Mayer /-bond for hard spheres equals —1 if 
the two spheres overlap and vanishes otherwise. For 
like species the deconvolution into scalar weight func- 
tions 2] can be written as fair) = — y^ 3 ,_ n w„ (R, r) * 

Here 



J2a'.T=o M ^' JPl UT n T (i,r), and given explicitly by 



U'3 



written as f u (r) = -^=o r- 



the spatial convolution of two functions is defined as 
(/ * .g)(x) = fdrf(r)g(x - r), and the metric M [l7| 
has components M^ v = M^ v = lifpt + i / = 3 and is 
zero otherwise. The crucial property of K(R,r) that we 
will exploit in the following is the group structure for the 
combined operation of matrix multiplication and (real- 
space) convolution. For the interactions between parti- 
cles of the same species K(2i?j,r) = K(Ri,r) * K(i?j,r), 
where the convolution product implies also matrix mul- 
tiplication, and the Mayer bond is just the special case 
f u (r) = -w 3 (2Ri,r) = -K 3 (2R l ,r). We exploit the 
fact that several matrices can be chained in or- 

der to model the interactions between unlike species via 
K(cr J ;j,r) = K(Ri,r) * K(dj,r) * K(dj,r) * K(i?j,r), where 
the lengthscales on the right hand side satisfy (p}. Fur- 
thermore the terms on the right hand side commute (i.e. 
the group is Abelian [16|). The Mayer bond between 
unlike species i and j is /y(r) = —K 3 (<7ij, 
hence be written as /y (r) = — , v T M w w 
K"{di,r) * K v T {dj,r) * w T (Rj,r). This identity, as well 
as that for the intra-species case above, can be verified 
by explicit algebra, most conveniently in Fourier space, 
where the convolutions become mere products and alge- 
braic theorems for trigonometric functions can be used 
to simplify the expressions. 



B. Algebraic-diagrammatic structure of the free 
energy functional 

Using the four weight functions w^(Ri,r), we build 
species-dependent weighted densities n M (i, r) in the stan- 
dard way via convolution with the bare density distribu- 
tion of the corresponding species, 



r) and can 
is(Ri,r)* 



n^(i,r) = w M (Ri,r) * pi{v). 



(2) 



We use the third-rank "junction" tensor J of Ref. [18[ 
in order to couple the n M (j, r) and hence generate terms 
that are non-linear in densities. Let us denote the com- 
ponents of J by J liVT . The tensor is symmetric under 
exchange of indices, and is non-zero only if /i + is + t = 6. 
One can specify J completely via the elements J 123 = 
J 033 = 1, and J 222 = l/(47r). As a basic building block 
for the construction of the density functional, we use 
the matrix N(i,r) of weighted densities 18j that is ob- 
tained by contracting the vector of weighted densities 
n T (i, r) with the J-tensor and lowering one of the indices 
via contraction with the metric M. Hence the compo- 
nents of the matrix N(i,r) are obtained as N^(i,r) = 



N(i,r) 



[n 3 (i,r) n 2 (i,r) m(i,r) n (z,r)^ 

n a (*,r) 23g£* ni(i>r ) 

n 3 (i,r) n 2 (i,r) 

V n a (i,r)/ 



(3) 



As an illustration of t he p ower of this formalized frame- 
work, we can obtain [18[ the additive FMT functional 
as the 03-component of J dx0od(X)i x ))i where the 
integration variable was renamed from r to x, and the 
zero-dimensional excess free energy is 



(rj) = {l-v) ln(l -V)+V=J2 



m(m — 1) 



> (4) 



with the (dummy) variable 7/ being the average occu- 
pation number of the zero-dimensional system. Here a 
function of a matrix is defined via its power series. In 
order to represent the mathematical structure of the in- 
tegrals in the density functional, we use the diagram- 
matic formulation of Ref. [Isj . see Fig. [2] for an overview. 
In particular, the star topology shown in Fig. con- 
stitutes the relevant type of diagram for additive hard 
sphere mixtures. The center of the diagram represents 
the integration variable x, the arms represent the weight 
functions w„(r), and the hllcd symbols represent the one- 
body density Pi(r) at space point(s) r. All spatial vari- 
ables are integrated over. The number of arms equals the 
order in density. Summing up all orders and using the 
coefficients of the power series of yields the Kierlik- 
Rosinberg free energy functional [18j. 

In the following we require only the last column of the 
Z-th power of the density matrix, obtained formally as 



^ ) ( ? ,r) = N'( ? , r) ■ (0,0,0, 1)* 



(5) 



where IM'(i, r) indicates the Z-th (matrix) power of IM(z, r), 
the dot represents the multiplication between a matrix 
and a vector and the superscript t indicates transposition. 
Explicitly, the four components in ([5]) are 



, I(*-l)(i-2) , 1-3 



ipl (i, r) = lnin l 3 



24tt 

-i Ki-i) 



"2 "3 



2„i-2 



8- 



-n 9 n 



2 "3 



(6) 

(7) 

(8) 
(9) 



where the superscripts of the weighted densities indi- 
cate (scalar) powers, the weighted densities are those of 
species i, i.e. n p = n^(i, r), and the arguments have been 
omitted for clarity. 

For non-additive mixtures we "transport" the expres- 
sions ©-([H]) via convolution with K(di,r). We hence 
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FIG. 2: Illustration of the topology of the diagrams that con- 
stitute the density series of the FMT free energy functional 
for various hard sphere mixtures. The one-body density dis- 
tribution of each species pi (r) (represented by a filled symbol) 
is connected by weight function bonds, w = (too, wi , W2, W3) 
(thin lines), to a central space integral (junction). The differ- 
ent topologies are for the binary additive hard spheres (a), bi- 
nary non-additive hard spheres (b,c) and ternary non-additive 
hard sphere mixtures (d). The numbers indicate the different 
species. In (b,c,d) the the junctions are joined by convolution 
kernel bonds K. The bonds that meet at a junction are mul- 
tiplied by third- rank tensors J. Multiplying with the metric 
M at the center of (c) joins the two K-bonds and restores the 
topology of (d). All junctions and all end points are integrated 



obtain a fourvector for each species i, indexed by /i = 
0, 1, 2, 3, and given as 

3 

*W(i,r) = £V(*,r)*4'>( i>r ). (1Q) 

These objects serve as ansatz functions for representing 
the free energy density; they are specific for each species 
(indicated by the argument i), are of Z-th order in density, 
and carry the dimension of (length) 3 "^. 



C. Rewriting the binary non-additive hard sphere 
functional 

Using the above definitions, we can write the den- 
sity functional for binary non-additive hard sphere mix- 
tures [ll| as 

F exc [p 1 ,p 2 ] = k B T [dy V M>^ ^ ( fc + '- 2 ) ! 



fX,v=Q 



k,l=0 



kin 



X$«(l,y)$(0 (2 , y ), 



(11) 



where we take the convention that the factorial vanishes 
for negative arguments and we have renamed the spa- 
tial integration variable from r to y. The scalar coef- 
ficients inside of the double sum over k, I are those in 



the Taylor expansion of the zero-dimensional excess free 
energy, which for a binary mixture is (^od^i + V2) — 
CVhiVi/Kk + l)(k + I - 1)]. Writing the coef- 
ficients as ( k + l )/[(k + l)(k + l - 1)] = (k + l- 2)1/ {kill) 
gives the form in (jlll) . 

A closed expression for the series (fTTj) can be obtained. 
This is identical to the previously given [ll[ form of the 
binary functional 

-Fexc[/9l,/92] = ksT 



fix I dx' K^{d 1 +d 2 ,x-x')^ liv {x,x!), 



(12) 



where the components of the free energy tensor are 
3 

$^(x,x')= \ a (l,x)A u p(2,x') 

a, 0=0 

x^od + ^^3(l,x)+n 3 (2,x')), (13) 

where ^odiv) = d a <J}oi(n)/dr] a is the a-th derivative of 
the Od excess free energy. Explicit expressions for the 
coefficients in ([TBI are 



n (i,x), 
ni(i,x)n 2 (i,x), 
[n 2 (i,x)] 3 
24tt ' 
ni(i,x), 

[n 2 (i,x)} 2 

n 2 (i,x), 
1, 



A 01 (i,x 
A 02 (i,x 

A 03 (i,x 
A n (i,x 

A 12 (i,x 
A 2 i(i,x 



(14) 



where i labels the species [i = 1, 2 for (fl~3|) ]. In the above 
we have exploited the convolution property 

K^{di+d 2 ,x-x') = 
3 

M TT , / dyK^{d u x-y)K T '»{d 2 ,x' -y). 

(15) 



3 

E 

t,t'=0 



Due to the group structure of the convolution kernels 
El , the identity K/ (d t + dj , r) = E'=o V (* , r ) * 
K T ' / (dj,r) holds, such that the lengthscale R\ 2 = o\ 2 — 
((T11 +o"22)/2 of the binary functional [111 is recovered as 
■R12 = di+d 2 . The structure of the diagrams correspond- 
ing to (fT2"|) is shown in Fig. and that corresponding to 
(ITT1) is shown in Fig. [5J:. 



D. Constructing a ternary non-additive hard 
sphere functional 

The benefit of re-writing the binary functional in the 
form (ITTT) is that this allows for straightforward general- 
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ization to three-component mixtures as 

3 

F cxc [pi,P2,P3] = k B T / rfx J ^ T 
x E ^ mi" *F (1 ' x >*® (2 ' (3 ' X) ' 

k,l,m— 

(16) 

I 



where 

$^ T (x,x',x") = 

3 

£ ^ IQ (l,x)^ /3 (2,x')A T7 (3,x") 

a, ,3,7=0 

x 4 Q d + ^ ] (n 3 (l,x) + n 3 (2,x') +n 3 (3,x")). (18) 

Eqs. ([HI) and ([17]) together with Eq. ([14]) (where i = 
1,2,3) prescribe the DFT for the general ternary hard 
sphere mixture. It is straightforward to verify that this 
functional reduces, in the corresponding limits, to the 
FMT for additive ternary har d sp here mixtures [l], 0, 
and for binary NAHS mixtures [11| . The diagrammatical 
structure of ([17]) is shown in Fig. [2jl. 



IV. RESULTS FOR BULK FLUID STRUCTURE 

We test the theory by calculating the partial two- 
body direct correlation functions for bulk fluids from 
Cii (|r-r'|) = -{k B T)^5^F cxc /5 Pl {v)8 P] {v')\ pk=const . 
Analytic expressions for the corresponding expressions 
Cij(q) in Fourier space can be obtained. Inserting these 
into the Ornstein-Zernike equation for ternary mixtures 
and Fourier transforming numerically yields partial pair 
correlation functions gij(r). Such results are shown in 
Fig. [3] for symmetric mixtures, an = 022 = 033 = 
with varying degree of non-additivity A12 = A23 = 
A23 = A, and at equal (and constant) bulk densities, 
pi = P2 = P3- We choose the overall packing fraction as 
77 = irpidfi/6 = 0.2 and compare to benchmark Monte 
Carlo simulation results for 1023 particles, obtained with 
10 6 attempted moves per particle of which the initial 10 5 
moves were used for equilibration. Except for (numeri- 
cally) small core violation, the OZ results reproduce the 
simulation data very well. This is also the case for the 



where (k + 1 + m — 2)!/(fc! l\ ml) is the coefficient of order 
ViVivT m the Taylor expansion of </>od(?7i +772+^3)- 
Again a closed expression can be obtained, which is in 
the form 



(17) 



slightly higher total packing fraction, rj — 0.25, where 
again MC and DFT data are compared in Fig. S] 
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1 1 1 1 1 1 1 

0.5 1 1.5 2 2.5 3 



r/a u 

FIG. 3: Partial pair correlation function gij(r) for symmetric 
ternary non-additive hard sphere mixtures with <tii = 022 = 
033 = a for varying degree of non-additivity A12 = A23 = 
A13 = A = —0.1,0,0.1 as a function of the scaled distance 
r/cr. Symbols denote Monte Carlo data; lines represent DFT 
results. Shown are the like correlation function (?n(r) (solid 
lines and crosses) and the unlike correlation function gi2(r) 
(dashed lines and pluses). The results for A = (0.1) are 
shifted upwards by two (four) units. 

In order to consider a fully asymmetric case, we have 
chosen the cry in order to mimick the correlations that 
were obtained in Molecular Dynamics simulations of an 
aqeuous binary salt solution [10| . Here the solvent is at 
large packing fraction and both types of ions are at van- 
ishing concentration. In this case the model parameters 
can be tuned to mimick the simulation data (not shown) . 
This constitutes a step towards possible applications of 
the current theory to electrolyte solutions. However, for 
mixtures with very asymmetric size ratios one would ex- 
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FIG. 4: Same as Fig. [3] but for total packing fraction r\ = 0.25. 

pect to find similar shortcomings as are present in the 
additive FMT [24| and indeed Percus-Yevick theory. 

V. CONCLUSIONS 

In conclusion, we have presented a fundamental mea- 
sure functional for ternary NAHS mixtures. The math- 
ematical structure of the functional is based on a dia- 
grammatic expansion in density with star-like (tree-like) 
topology of the diagrams. We have shown that pair corre- 
lation functions obtained from the density functional via 
the Ornstein-Zernike route agree well with computer sim- 



ulation data. Possible applications of the current theory 
include its use as a reference system in modeling elec- 
trolytes in bulk and in inhomogeneous situations. Fur- 
thermore, considering one-dimensional cases, see Ref. [25[ 
for the exact solution of bulk properties of binary mix- 
tures, could be interesting. 

It is worthwhile to discuss possible generalization to 
mixtures with more than three components. The gen- 
eral ternary case presented here rests crucially on the 
decomposition of the interaction distances into suit- 
able length scales Ri and di, cf. ([1]). This allowed to 
use suitable diagrams with tree-like structure, cf. Fig.[2}i. 
However, (p} does not have a general solution (for Ri and 
di once the cry are prescribed) for M > 4 components. 
That this is true can be gleaned from the fact that the cr^ 
are M(M + l)/2 independent constants, whereas the Ri 
and di constitute only 2M free parameters. These num- 
bers, seemingly by accident, match in the special case 
AI = 3. Hence we can conclude that possible extensions 
to four and more components requires further changes in 
the mathematical structure of the present FMT. 
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